Computer simulations of history of life: speciation, emergence of complex species 

from simpler organisms, and extinctions 
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We propose a generic model of eco-systems, with a hierarchical food web structure. In our 
computer simulations we let the eco-system evolve continuously for so long that that we can monitor 
extinctions as well as speciations over geological time scales. Speciation leads not only to horizontal 
diversification of species at any given trophic level but also to vertical bio-diversity that accounts for 
the emergence of complex species from simpler forms of life. We find that five or six trophic levels 
appear as the eco-system evolves for sufficiently long time, starting initially from just one single 
level. Moreover, the time intervals between the successive collections of ecological data is so short 
that we could also study "micro" -evolution of the eco-system, i.e., the birth, ageing and death of 
individual organisms. 



I. INTRODUCTION 

During last one decade, theoretical research on co- 
evolution of species in eco-systems and the statistics of 
extinctions have been strongly infleunced by the pioneer- 
ing interdisciplinary works of Per Bak and his collabo- 
rators PiSSH- In the same spirit, we address some 
fundamental questions of evolutionary ecology from the 
perspective of statistical physics. 

How did the higher species emerge in eco-systems in- 
habited initially only by primitive forms of life, like bac- 
teria and plankton 0,|a,l3|'^ The available record of the 
history of life, written on stone in the form of fossils, 
is incomplete and ambiguous. An alternative enterprise 
seeks to recreate the evolution on a computer by simu- 
lating theoretical models 0, 13 • In this paper we propose 
a theoretical model that not only addresses the ques- 
tion raised above but also provides a versatile concep- 
tual tool for studying evolutionary ecology. In particu- 
lar, it describes both the "macro" -evolutionary processes 
(e.g., origin, evolution and extinction of species) as well 
as "micro" -evolution, (e.g., age-distribution in the popu- 
lation of a species, mortality rates, etc.). 

If watched over a short period of time, the dynam- 
ics of the eco-system appears to be dominated by birth 
and death of the individual organisms as well as by the 
prey-predator interactions. However, over longer pe- 
riod of time, one would see not only extinction of some 
species but also the appearance of new ones. Besides, 
in many situations, macro-evolutionary changes occur at 
rates that are comparable to those of the ecological pro- 
cesses 0, ^1 . The artificial separation of this process 
into "ecological" time scales and "geological" time scales 
|l2j | has been made in many earlier theoretical works only 
for the convenience of modelling. 

The "ecological" models, that describe population 
dynamics in detail using, for example, the Lotka- 
Volterra equations [TsLIT^I usually ignore the slow macro- 
evolutionary changes in the eco-system; hardly any ef- 
fects of these would be observable before the computer 
simulations would run out of computer time [l5l |. On 



the other hand, in order to simulate the billion-year old 
history of life on earth with a computer, the elementary 
time steps in "evolutionary" models have to correspond 
to thousands of years, if not millions; consequently, the 
finer details of the ecological processes over shorter pe- 
riods of time cannot be accounted for by these models 
in any explicit manner Ill.|l6|. Limitations of these ap- 
proaches are well knowriflTl. Il8l ITgj . Moreover, most of 
the recent computer models of ageing focus attention 
on only one isolated species and, therefore, cannot cap- 
ture macro-evolutionary phenomena like, for example, 
extinctions which depend crucially on the prey-predator 
interactions. 

We wish to develop one single theoretical model which 
would be able to describe the entire dynamics of an eco- 
system since the first appearance of life in it up till now 
and in as much detail as possible. This dream has now 
come closer to reality, ma inly because of the availability 
of fast computers [H IH IM El IM El ■ It has become 
feasible now to carry out computer simulations (in-silico 
experiments) of eco-system models where, each time step 
would correspond to typical times for "micro" -evolution 
while each of the simulations is run long enough to cap- 
ture "macro" -evolution. 

The prey-predator relations in any eco-system are 
usually described graphically in terms of food webs 
(271,1^ mill. More precisely, a food web is a directed 
graph where each node is labelled by a species' name 
and each directed link indicates the direction of flow of 
nutrient (i.e., from a prey to one of its predators). We 
incorporate in our model the hierarchical organization of 
the species at different trophic levels of the food web. 

In real eco-systems, the food web is a slowly evolv- 
ing dynamic network. For example, species are known 
to change their food habits "si'l. These changes in diets 
may be caused by scarcity of the normal food and abun- 
dance of alternative food resources. Moreover, higher 
organisms appear through speciation in an eco-system 
that initially began with only simple forms of life. These 
not only occupy new trophic levels but also introduce 
new prey-predator interactions with the existing species. 
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Therefore, it is also desirable that these self-organizing 
features of natural eco-systems should be reproduced, at 
least qualitatively, by the theoretical models. 

The aim of this paper is to propose a model that would 
capture the desirable features of eco-systems outlined 
above. Higgs, McKane and collaborators 0, 0| have 
developed a model, called the Webworld model, which 
was aimed at linking the ecological modeling of food web 
architecture with the evolutionary modeling of speciation 
and extinction. The spirit of our model is very similar 
although the details of the mathematical formulation of 
the two models are quite different. 

II. THE MODEL 

We model the eco-sytem as a dynamic hierarchical 
network. The "micro" -evolution, i.e., the birth, growth 
(ageing) and natural death of the individual organisms, 
in our model is captured by the intra-node dynamics. 
The "macro "-evolution, e.g., adaptive co-evolution of the 
species, is incorporated in the same model through a 
slower evolution of the network itself over longer time 
scales. Moreover, as the model eco-system evolves with 
time, extinction of species is indicated by vanishing of the 
corresponding population; thus, the number of species 
and the trophic levels in the model eco-system can fluc- 
tuate with time. Furthermore, the natural process of 
speciation is implemented by allowing re-occupation of 
vacant nodes by mutated versions of non-extinct species. 

A. Architecture of the network 

Each node of the network represents a niche that can 
be occupied by at most one species at a time. The num- 
ber of nodes in the trophic level £ is where m is 
a positive integer. We assume only one single species 
at the highest level £ = 1. The allowed range of £ is 
1 < £ < £max(t), where £max{t) is a time-dependent num- 
ber in our new model. In other words, in contrast to all 
cited earlier models, the numerical value of £max in our 
new model is not put in by hand, but is an emergent 
property of the eco-system. 

B. Prey- predator interactions and intra-species 
competitions 

The prey-predator interaction between two species 
that occupy the nodes i and k at two adjacent trophic 
levels is represented by Jik\ the three possible values of 
Jik are ±1 and 0. The sign of Jik indicates the direction 
of trophic flow, i.e. from the lower to the higher level. 
Jik is -1-1 if I eats k and it is —1 if A: eats i. If there is 
no prey-predator relation between the two species i and 
fc, we must have Jik = 0. Although there is no direct 
interaction between species at the same trophic level in 



our model, they can compete, albeit indirectly, with each 
other for the same food resources available in the form of 
prey at the next lower trophic level. 

We now argue that the elements of the matrix J ac- 
count not only for the mter-species interactions but also 
for the mfra-species interactions arising from the compe- 
tition of individual organisms for the same food resources. 
Let be the number of all prey individuals for species 
i on the lower trophic level, and S~ be m times the num- 
ber of all predator individuals on the higher trophic level. 
Since we assume that a predator eats m prey per time 
interval, Sf gives the amount of total food available for 
species i, and S~ is the total contribution of species i 
to the pool of food required for all the predators on the 
higher level. If the available food Sf is less than the 
requirement, then some organisms of the species i will 
die of starvation^ even if none of them is killed by any 
predator. 

The mira-species competition among the organisms 
of the same species for limited availability of resources, 
other than food, imposes an upper limit Umax of 
the allowed population of each species; Umax is time- 
independent parameter in the model. Thus, the total 
number of organisms nit) at time t is given by n{t) = 

Erjs.it). 

If Hi — is larger than then food shortage will 
be the dominant cause of premature death of a fraction 
of the existing population of the species i. On the other 
hand, if S~ > Ui — , then a fraction of the existing 
population will be wiped out primarily by the predators. 

In order to capture the starvation deaths and killing 
by the predators, in addition to the natural death due to 
ageing, a reduction of the population by 

C max{S^ , Ui — S^) (1) 

is implemented at every time step, where Ui is the pop- 
ulation of the species i that survives after the natural 
death. C is a constant of proportionality. If this leads to 
Ui < 0, species i becomes extinct. 

We assume that the simplest species occupying the 
lowest trophic level always get enough resources that nei- 
ther natural death nor predators can affect their popula- 
tion. 



C. Collective characteristics of species 

An arbitrary species i is collectively characterized by 

M 

(i) the minimum reproduction age X^^p (z) , 

(ii) the birth rate M{i), 

(iii) the maximum possible age Xmaxi^) that depends 
only on the trophic level occupied by the species. 

An individual of the species i can reproduce only after 
attaining the age Xj-epii)- Whenever an organism of this 
species gives birth to offsprings, M (i) of these are born si- 
multaneously. None of the individuals of this species can 
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live longer than Xmax{i)i even if an individual manages 
to escape its predators or starvation. 

D. Mutations 

With probability Pmut per unit time, each of the 
species randomly increases or decreases, with equal prob- 
ability, their Xrep and M by unity. {X^ep is restricted 
to remain in the interval from 1 to Xmax, and M > 0.) 
Moreover, with the same probability Pmut per unit time, 
they also re-adjust one of the links J from prey and one 
of the links J to predators [s^l- If the link Jij to the 
species i from a higher level species j is non-zero, it is as- 
signed a new value of Jij = Jji = 0. On the other hand, 
if the link Jik to a species i from a lower species k is zero, 
the new values assigned are Jik = l,Jki = ^1- These re- 
adjustments of the incoming and outgoing (in the sense 
of nutrient flow) interactions are intended to capture the 
facts that each species tries to minimize predators but 
look for new food resources. 



E. Speciation 

The niches (nodes) left empty because of extinction are 
re- filled by new species, with probability per unit time 
■ AH the simultaneously re- filled nodes in a trophic 
level of the network originate from one common ancestor 
which is picked up randomly from among the non-extinct 
species at the same trophic level. All the interactions J 
of the new species are identical to those of their com- 
mon ancestor. The characteristic parameters Xrep, M of 
each of the new species differ randomly by ±1 from the 
corresponding parameters for their ancestor. 

However, occasionally, all the niches at a level may lie 
vacant. Under such circumstances, all these vacant nodes 
are to be filled by a mutant of the non-extinct species 
occupying the closest lower populated level. As stated 
above, the lowest level, that is populated by the simplest 
species, never goes extinct; the possible ageing of the 
species at the lowest level |3J| is not relevant here. All the 
individual organisms of the new species are assumed to 
be newborn babies that begin ageing with time just like 
the other species. Since space does not enter explicitly 
in our model, it does not distinguish between sympatric 
and allopatric speciation [? ]. 

F. Emergence of new trophic level 

In order to understand why the total number of trophic 
levels in food webs usually lie between 4 and 6, we al- 
lowed adding a new trophic level to the food web, with 
a small probability piev per unit time, provided the to- 
tal bio-mass distributed over all the levels (including the 
new one) does not exceed the total bio-mass available in 
the eco-system. This step is motivated by the fact that 



real ecosystems can exhibit growing bio-diversity over 
sufficiently long period of time. Increase of the number 
trophic level means the diversification at the erstwhile 
topmost level as well as all the lower levels and the emer- 
gence of yet another dominating species that occupies 
the new highest level. The total number of levels, which 
determines the lengths of the food chains, depends on 
several factors, including the available bio-mass 36]. 

G. Birth and natural death of organisms 

At each time step, each individual organism a of the 
species i gives birth asexually to M{i) offsprings with a 
probability pb(*, a)- We also assume the time- dependent 
probability pb(i,a) is a product of two factors. One of 
these two factors decreases linearly with age, from unity, 
attainable at the minimum reproduction age, to zero at 
the maximum lifespan. The other factor is a standard 
Verhulst factor 1 — ni/n„iax which takes into account the 
fact that the eco-system can support only a maximum of 
n-max individual organisms of each species. Thus, Pbii, ct) 
is equal to the Verhulst factor at X = Xrep- 

Each individual organism, irrespective of its age, can 
meet its natural death. However, the probability pd of 
this natural death depends on the age of the individ- 
ual. In order to mimic age-independent constant mor- 
tality rate in childhood, we assume the probability pd 
of "natural" death (due to ageing) to be a constant 
Pd — &yi^\—r{Xmax — Xrep)/M]^ (with a small fraction 
r), so long as X < Xrep- However, for X > Xrep, the 
probability of natural death is assumed to increase fol- 
lowing the Gompertz law pd — exp[—r{Xrnax — X)/M]. 
Note that, for a given Xmax and Xrep, the larger is the M 
the higher is the pd for any age X. Therefore, in order 
maximize reproductive success, each species has a ten- 
dency to increase M for giving birth to larger number of 
offsprings whereas the higher mortality for higher M op- 
poses this tendency js^l- However, even with a constant 
Pd = 0.1 wc found qualitatively similar results. 

H. Summary of the dynamics of the eco-system 

The state of the system is updated in discrete time 
steps where each step consists of a sequence of six stages: 

/- Birth 

II- Natural death 

III- Mutation 

IV- Starvation death and killing by prey 

V- Speciation 

VI- Emergence of new trophic level 

In all our simulations we began with random initial 
condition, except for M = 1 for all species, mostly with 
only three levels in the food web, and let the eco-system 
evolve for T^, time steps before we started collecting eco- 
logical and evolutionary data from it; these data were 
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FIG. 1: Log-log plots of the distributions of the lifetimes 
of the species in an eco-system with the total biomass of 
4nmaa;. The symbols +, x, * and □ correspond, respec- 
tively, to 10*, 10^, 10'' and 10* time steps. Each of the 
data points have been obtained by averaging over 6400 sam- 
ples 640 samples (x), 64 samples (*) and 1 sample (□). 
The line with slope —2 corresponds to a power law distribu- 
tion that has been predicted by many theories. The com- 
mon parameters for all the plots are m = 2, Umax = 100, 
Psp = 0.1,p^„t = 0.001, = 0.0001, C = 0.2, r = 0.05. 
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FIG. 2: Semi-log plots of the distributions of the minimum 
reproductive age Xre.p of the species. The symbols x, *, □ 
and ■ correspond, respectively, to 10*, 10^, 10®, 10^ and 10* 
time steps. Each of the data points have been obtained by 
averaging over 6400 samples 640 samples (x), 64 samples 
(*) and 1 sample (□ and ■). The values of the common 
parameters for all the plots are identical to those used in figQ 



collected for the subsequent T = 5Tm time steps where 
the longest runs were for T = 10^. We have not observed 
any quahtative differences in the data for rimax = 100 
and rimax — 1000, keeping all the other parameters 
same. Most of our simulations were carried out with 
m = 2, as we did not observe quahtative differences be- 
tween the data for m — 2 and m — 3,4 in test runs. 
The maximum lifespans in the levels were assumed to be 
Xmax = 100,71,50,35,25,.... starting from the highest 
level. 



III. RESULTS 

A. Lifetime distributions 

Several theories, based on extremely simple models, 
claim that the distribution of the lifetimes of the species 
should follow a power law with a slope of —2 on the log- 
log plot. The distributions of the lifetimes of the species 
in our model are shown in fig^for a few different sets of 
parameter values. Although our data do not rule out an 
approximate power law over limited regime of lifetimes, 
one single power law over the entire range of lifetimes 
seems impossible. 



B. Distribution of minimum reproductive ages 

In fig[21 we plot the distributions of the minimum re- 
productive age Xrep of the species for several different 
sets of values of the model parameters. Although over 



relatively short time scales of observation this distribu- 
tion appears quite broad it narrows down with evolution 
and the non-zero values of this distribution correspond 
to reasonable values of age. 



C. The number of trophic levels 

Due to the randomness in the evolutionary process, 
occassionally, all of the niches in a level (except the low- 
est one) may lie vacant. We have monitored imax(t) and 
also A/'(t), the number of those levels at time t in which 
at least one niche is occupied by a non-extinct species. 
In figOl we plot A/" as a function of time for one single 
run. This clearly shows how, over geological time scales, 
(^max reaches 6. In this run, the sixth level (the highest 
one) emerges after 10^ time steps. It also demonstrates 
that at all stages of evolution, the number J\f{t) keeps 
fluctuating. During the very late stages, Af keeps fluctu- 
ating between 5 and 6, although Af is more often 6 than 
5 for all times beyond 10^. The ratio of occurrences of 
six levels and five levels in the eco-system stabilized only 
after 10^ time steps. 

We have also computed the distributions (histograms) 
of J\f by averaging the data over large number of runs. 
As shown in fig01 the distribution becomes narrower for 
longer runs and the trend indicates than in the extreme 
long time limit would be sharply peaked around one sin- 
gle value of Af, as indicated by the figOl 
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FIG. 3: Semi-log plot of the number M{t) of the trophic levels 
with at least one non-extinct species plotted against time t in 
one single run. The values of the parameters are identical to 
those used in fig0 
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FIG. 4: Distribution (histogram) of the number J\f of the 
trophic levels with at least one non-extinct species. Identical 
symbols in this figure and in figQ correspond to the same 
runs. The data for the longest time steps are connected by 
lines to emphasize the shrinking trend of the distribution with 
time. 



IV. SUMMARY AND CONCLUSION 

In this paper we have introduced a theoretical model 
of eco-systems with a generic hierarchical trophic level 
structure. Because of data collection at sufFciently short 
intervals, we have been able to monitor the ecological 
phenomena like, for example, birth, ageing and death of 
individual organisms and, hence, the population dynam- 
ics of the species. We have also been able to run our 
simulations upto sufficiently long times (10^, with sta- 
tionarity achieved at around 10^) so that the model also 



accounted for macro-evolutionary phenomena like extinc- 
tions of species as well as speciation that leads not only 
to emergence of species at the existing levels of complex- 
ity but also to higher species that occupy an altogether 
new trophic level in the food web. 

From the infinite possible life forms, we start with one 
or a few, and then let our ecosystem grow in diversity and 
complexity until the limitations of biomass restrict it to 
hundreds of individuals in dozens of species, organized 
into about five trophic levels. 

Although our model food web is hierarchical, it is not 
a tree-like structure. The hierarchical architecture helps 
us in capturing a well known fact that in the normal 
ecosystems the higher is the trophic level the fewer are 
the number of species. It is well known that the body 
size and abundance of a species are strongly correlated 
to their positions as well as to their interactions with 
other species in the food web [sM l39l | . If we neglect par- 
asites and herbivorous insects on trees, then, in general, 
predators are fewer in number and bigger in size as com- 
pared to their prey species js^. This is very naturally 
incorporated in the hierarchical food web structure of our 
model. Let us assume that in the model the body size 
of individual organisms on each level £ is about m times 
smaller than that on its predator level ^ — 1. On the other 
hand, the maximum possible populations of organisms, 
including all the nodes, in a level € is to times that at 
the level ^ — 1. Consequently, the maximum amount of 
biomass on each level is, approximately, the same. 

Since each individual organism appears explicitly in 
our model, one could, at least in principle, assign a 
genome to each individual and describe Darwinian selec- 
tion which takes place at the level of organisms |3 [43| . 
Unfortunately, additional ad-hoc assumptions would be 
required to relate the genome with the reproductive suc- 
cess I^H^. Instead of introducing an ad- hoc mathemat- 
ical formula to relate genotype with phenotype, we have 
worked directly with phenotype, particularly, quantities 
that decide the reproductive success of the organisms; 
these quantities are X^ep, Xmax and M. 

From the perspective of self-organization, the new 
model surpasses all cited previous models as not only 
the characteristic collective properties of the species but 
even the nature of inter-species interactions as well as 
the total number of trophic levels in the food web are 
determined by self-organization of the eco-system. 
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